%pylab inline
colors = [[0.1,0.6,0.5], [0.9,0.2,0.5], [0.9,0.5,0.1]]
Populating the interactive namespace from numpy and matplotlib
from PyDynamic.uncertainty.propagate_DFT import GUM_DFT,GUM_iDFT
from PyDynamic.uncertainty.propagate_DFT import DFT_deconv, AmpPhase2DFT
from PyDynamic.uncertainty.propagate_DFT import DFT_multiply
#%% reference data
ref_file = np.loadtxt("DFTdeconv reference_signal.dat")
time = ref_file[:,0]
ref_data = ref_file[:,1]
Ts = 2e-9
N = len(time)
#%% hydrophone calibration data
calib = np.loadtxt("DFTdeconv calibration.dat")
f = calib[:,0]
FR = calib[:,1]*np.exp(1j*calib[:,3])
Nf = 2*(len(f)-1)
uAmp = calib[:,2]
uPhas= calib[:,4]
UAP = np.r_[uAmp,uPhas*np.pi/180]**2
#%% measured hydrophone output signal
meas = np.loadtxt("DFTdeconv measured_signal.dat")
y = meas[:,1]
# assumed noise std
noise_std = 4e-4
Uy = noise_std**2
Consider knowledge about the measurement system is available in terms of its frequency response with uncertainties associated with amplitude and phase values.
H = (|H(f1)|, …, ∠H(fN))
uH = (u|H(f1)|, …, u∠H(fN))
figure(figsize=(16,8))
errorbar(f * 1e-6, abs(FR), 2 * sqrt(UAP[:len(UAP) // 2]), fmt= ".-", alpha=0.2, color=colors[0])
xlim(0.5, 80)
ylim(0.04, 0.24)
xlabel("frequency / MHz", fontsize=22); tick_params(which= "both", labelsize=18)
ylabel("amplitude / V/MPa", fontsize=22);
figure(figsize=(16,8))
errorbar(f * 1e-6, unwrap(angle(FR)) * pi / 180, 2 * UAP[len(UAP) // 2:], fmt= ".-", alpha=0.2, color=colors[0])
xlim(0.5, 80)
ylim(-0.2, 0.3)
xlabel("frequency / MHz", fontsize=22); tick_params(which= "both", labelsize=18)
ylim(-1,1)
ylabel("phase / rad", fontsize=22);
The measurand is the input signal x = (x1, …, xM) to the measurement system with corresponding measurement model given by
y[n] = (h * x)[n] + ε[n]
Input estimation is here to be considered in the Fourier domain.
The estimation model equation is thus given by
with - Y(f) the DFT of the measured system output signal -HL(f) the frequency response of a low-pass filter
Estimation steps
- DFT of y and propagation of uncertainties to the frequency domain
- Propagation of uncertainties associated with amplitude and phase of system to corr. real and imaginary parts
- Division in the frequency domain and propagation of uncertainties
- Multiplication with low-pass filter and propagation of uncertainties
- Inverse DFT and propagation of uncertainties to the time domain
With the DFT defined as
with βn = 2πn/N, the uncertainty associated with the DFT outcome represented in terms of real and imaginary parts, is given by
Y,UY = GUM_DFT(y,Uy,N=Nf)
figure(figsize=(18,6))
subplot(121)
errorbar(time*1e6, y, sqrt(Uy)*ones_like(y),fmt=".-")
xlabel("time / µs",fontsize=20); ylabel("pressure / Bar",fontsize=20)
subplot(122)
errorbar(f*1e-6, Y[:len(f)],sqrt(UY[:len(f)]),label="real part")
errorbar(f*1e-6, Y[len(f):],sqrt(UY[len(f):]),label="imaginary part")
legend()
xlabel("frequency / MHz",fontsize=20); ylabel("amplitude / au",fontsize=20);
In practice, the frequency response of the measurement system is characterised in terms of its amplitude and phase values at a certain set of frequencies. GUM uncertainty evaluation, however, requires a representation by real and imaginary parts.
Hk = Akcos (Pk) + jAksin (Pk)
GUM uncertainty propagation
H, UH = AmpPhase2DFT(np.abs(FR),np.angle(FR),UAP)
Nf = len(f)
figure(figsize=(18,6))
subplot(121)
errorbar(f*1e-6, H[:Nf], sqrt(diag(UH[:Nf,:Nf])),fmt=".-",color=colors[2],alpha=0.2)
xlabel("frequency / MHz",fontsize=20); ylabel("real part / au",fontsize=20)
subplot(122)
errorbar(f*1e-6, H[Nf:],sqrt(diag(UH[Nf:,Nf:])),fmt=".-",color=colors[2],alpha=0.2)
xlabel("frequency / MHz",fontsize=20); ylabel("imaginary part / au",fontsize=20);
The deconvolution problem can be decomposed into a division by the system's frequency response followed by a multiplication by a low-pass filter frequency response.
which in real and imaginary part becomes
Sensitivities for division part
Uncertainty blocks for multiplication part
# low-pass filter for deconvolution
def lowpass(f,fcut=80e6):
return 1/(1+1j*f/fcut)**2
HLc = lowpass(f)
HL = np.r_[np.real(HLc), np.imag(HLc)]
XH,UXH = DFT_deconv(H,Y,UH,UY)
XH, UXH = DFT_multiply(XH, UXH, HL)
figure(figsize=(18,6))
subplot(121)
errorbar(f*1e-6, XH[:Nf], sqrt(diag(UXH[:Nf,:Nf])),fmt=".-",color=colors[2],alpha=0.2)
xlabel("frequency / MHz",fontsize=20); ylabel("real part / au",fontsize=20)
subplot(122)
errorbar(f*1e-6, XH[Nf:],sqrt(diag(UXH[Nf:,Nf:])),fmt=".-",color=colors[2],alpha=0.2)
xlabel("frequency / MHz",fontsize=20); ylabel("imaginary part / au",fontsize=20);
The inverse DFT equation is given by
The sensitivities for the GUM propagation of uncertainties are then
GUM uncertainty propagation for the inverse DFT
xh,Uxh = GUM_iDFT(XH,UXH,Nx=N)
ux = np.sqrt(np.diag(Uxh))
figure(figsize=(16,8))
plot(time*1e6,xh,label="estimated pressure signal",linewidth=2,color=colors[0])
plot(time * 1e6, ref_data, "--", label= "reference data", linewidth=2, color=colors[1])
fill_between(time * 1e6, xh + 2 * ux, xh - 2 * ux, alpha=0.2, color=colors[0])
xlabel("time / µs", fontsize=22)
ylabel("signal amplitude / MPa", fontsize=22)
tick_params(which= "major", labelsize=18)
legend(loc= "upper left", fontsize=18, fancybox=True)
xlim(0, 2);
figure(figsize=(16,8))
plot(time * 1e6, ux, label= "uncertainty", linewidth=2, color=colors[0])
xlabel("time / µs", fontsize=22)
ylabel("uncertainty / MPa", fontsize=22)
tick_params(which= "major", labelsize=18)
xlim(0,2);
Y,UY = GUM_DFT(y,Uy,N=Nf)
H, UH = AmpPhase2DFT(A, P, UAP)
XH,UXH = DFT_deconv(H,Y,UH,UY)
XH, UXH = DFT_multiply(XH, UXH, HL)